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We introduce a system of two linearly coupled discrete nonlinear Schrodinger equations (DNLSEs) , 
with the coupling constant subject to a rapid temporal modulation. The model can be realized in 
bimodal Bose-Einstein condensates (BEC). Using an averaging procedure based on the multiscale 
method, we derive a system of averaged (autonomous) equations, which take the form of coupled 
DNLSEs with additional nonlinear coupling terms of the four-wave-mixing type. We identify stabil- 
^ ' ity regions for fundamental onsite discrete symmetric solitons (single-site modes with equal norms 

O ' in both components), as well as for two-site in-phase and twisted modes, the in-phase ones be- 

, ing completely unstable. The symmetry-breaking bifurcation, which destabilizes the fundamental 

symmetric solitons and gives rise to their asymmetric counterparts, is investigated too. It is demon- 
' strated that the averaged equations provide a good approximation in all the cases. In particular, 

the symmetry-breaking bifurcation, which is of the pitchfork type in the framework of the averaged 
equations, corresponds to a Hopf bifurcation in terms of the original system. 
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PACS numbers: 42.65.-k, 42.50.Ar, 42. 81. Dp 



I. INTRODUCTION 



The problem of the influence of rapidly varying perturbations on solitons belongs to the general topic of the "soliton 
management" Along with other problems of this kind, this one has drawn considerable attention. The motivation 
is to investigate new nontrivial dynamics that may be induced by rapidly varying perturbations, and predict, in 
this way, new types of solitons in such systems. In terms of mechanical systems, with few degrees of freedom, the 
\^ closest counterparts of this setting are presented by the Kapitza pendulum and the motion of a charged particle in 
rapidly oscillating electromagnetic fields. In either case, the averaged dynamics is governed by an effective potential, 
' \ different from the original one, whose fixed points may give rise to novel dynamical states. For the discrete nonlinear 
I . Schrodinger (DNLS) equation, problems of this kind ("rapid management") were considered for a variable second- 
' order discrete dispersion jS] , and for the variable nonlinearity [1] . The respective averaged equations take the form of 
T— I I a nonlocal generalization of the DNLSE, which gives rise to new types of localized breathers. Recently, the existence 
' of a soliton in the DNLSE with an external drive has been demonstrated in Ref . [5| , and the influence of a temporal 
. ^ ' delay in the onsite nonlinearity on the self-trapping of discrete solitons was studied in Ref. @ . 

• The above-mentioned examples pertain to single-component discrete systems. A natural physically relevant gen- 
J-j \ eralization is to explore the influence of rapidly varying strong perturbations on the dynamics of discrete vectorial 
. . . ■ (two-component) solitons. In this work,we introduce a system of two linearly coupled DNLSEs with a rapidly varying 
coupling parameter. This example, which is interesting in its own right, is a straightforward model of Bose-Einstein 
condensates (BECs) confined in two parallel tunnel-coupled cigar-shaped traps, combined with a deep optical lattice 
(OL), in a case when the linear-coupling parameter is subject to the temporal modulation 0, Another imple- 
mentation of the model is possible in terms of the binary BEC (also loaded into a deep OL) with two components 
linearly coupled by a resonant electromagnetic wave [l^ , whose amplitude may also be periodically modulated in 
time In the latter case, the coupled system of the Gross-Pitaevskii equations also contains nonlinear-interaction 
terms, accounting for collisions between atoms belonging to the different species. However, the coefficient in front of 
the latter terms may be effectively switched off by means of the Feshbach resonance affecting the inter-species colli- 
sions (see, e.g., Ref. [HI), which we assume below. Our goal is to derive effective averaged equations approximating 
this model, and thus predict discrete solitons in it~both symmetric ones, and asymmetric states generated by the 
symmetry-breaking bifurcation in the dual-core system. The analysis of this bifurcation in "unmanaged" dual-core 
DNLSE systems, with a constant linear-coupling constant, was developed in Refs. [l^ and (Tsj . 

The paper is structured as follows. The model is introduced in Section II, where we also derive the approximation 
based on the averaged equations. Symmetric discrete solitons, of both single-site and two-site types (in-phase and 
twisted ones, in the latter case), are investigated in Section III. In Section IV, we explore the symmetry-breaking bifur- 
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cation, which destabihzes the single-site symmetric sohtons, giving rise to their asymmetric counterparts. Conclusions 
are formulated in Section V. 



II. THE MODEL AND AVERAGED EQUATIONS 



Following the above discussion, we introduce the following system of linearly coupled DNLSEs, written in a scaled 
form, 

«^Wn + Kl(w„+1 - M„-l) +7o|w„|^W„ + /(O^n = 0, (la) 

at 

i-r:Vn+ H2{Vn+l+Vn-l)+la\Vn?Vn + I{t)Un = 0, (lb) 

at 

for amplitudes w„(i) and w„(i) of the two tunnel-coupled BECs trapped in the deep OL [l2| . with time-dependent 
coupling coefficient f{t). Here, we consider the case of the strong high-frequency temporal modulation, f{t) — 
fo + fi sin(r2i),with /i ~ > 1- 

Different methods can be used to derive averaged equations for slowly varying fields, such as the multi-scale 
expansion |il4;] and the Kapitza method [2]. In the following, we derive the averaged system by way of the former 

technique. To this end, we first define C(r) = /i sin(r), where t = fit and /i = /i/fi, so that /i = 0{1). To remove 
the secular part of the coupling terms, we use a linear transformation. 



\ / cos(C_i(t)) -isin(C_i(T)) 
Vn J I -isin(C_i(T)) cos(C_i(t)) J\lpn 



(2) 



where, in the case of the general time modulation, C_i(t) = C{t') dr' — P ^ C{t') dr' dr, and P is the 

period of C(t). For the periodic coupling defined above, C_i = — /i cos{ilt), and P = 2tt. 

Substituting (/)„ = $„(t) -I- 0{l/fl) and 4'n = ^n(i) + 0(1/^^) and averaging over the fast time variable, r, it is 
straightforward to derive the following averaged equation for $„ and , 

z^$„ + Ki($„+i + $„-i) + ? [(3 + a)|$„P + 2(1 - a)\^n\^] + ^(a - + fo^n = 0, (3a) 

dt 4 4 

+ ^^{'fn+i + *„-i) + ^ [2(1 - a)|$„P + (3 + a)|«'„n + ^{a - + /o$„ = 0, (3b) 

dt 4 4 

where a — P^^ cos(4C_i(t)) dr. For the form of the periodic modulation adopted above, 

<j={2n)-^J^ cos (4/1 cos(r)) dr = Jo (4/1) , (4) 

where Jq is the Bessel function. Note that the nonlinear coupling terms generated by the averaging method in Eqs. 
are sensitive to the phase difference between the two fields, i.e., these terms are of the four- wave-mixing type. 



III. SYMMETRIC LOCALIZED MODES 



A. The general approach 

To check the accuracy of our averaging analysis, in the following we explore the existence and stability of discrete 
solitons predicted by Eqs. ([T]) and ([3]). First, we consider stationary solutions which are symmetric with respect to 
the coupled subsystems, i.e., with — "^n- 

Using the Newton- Raphson continuation method, we look for stationary solutions to Eqs. In this case, the 
solution is sought for in the form of $„(t) = e^*t/„ and 4'„(t) = e^^Vm where the frequency is assumed to be fixed as 
A = 1, by means of an obvious rescaling. Substituting this in Eqs. ([3]), we arrive at the following equations for the 
real- valued functions, [/„ and Vn'- 

-A[/„-HKi([/„+i + C/„_i) + i7o((3 + fT);72 + (l-a)V;2)[/„ + /oi4 = 0, (5a) 

-KV^ + K2{Vn+l+Vn-l) + \la{{l~a)Ul + {i + a)V^)Vn + hUn = 0. (5b) 
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To test the stability of the stationary sohitions, we substitute the usual expression for a weakly perturbed solution, 
^n{t) = e"^' (Un + e*^*M„) and — e'^* (Vn + e*^'w„), into Eqs. ^ and perform the linearization of the equations, 

arriving at the eigenvalue problem for perturbation frequencies A, 

Ki(u„+i + + ^ [(3 + cT)U^i2un + K) + (1 - ^) {V^i'^Un " <) + 2U„VnK)] + /o"n = (A + A) u„, (6a) 

At2(«„+i+S„-i) + ^ [(3 + ^)K?(2«n + «;j + (l-fT)(t/2(2f}„-S;) + 2t/„K<)] +/o^^ - (A + A)S„. (6b) 

The discrete soliton is stable if all eigenvalues A have Im(A) = 0. 



B. Fundamental discrete solitons (single-site modes) 

We start with the consideration of the existence and stability of a symmetric single-site mode, which, in the 
anticontinuum limit of ki — K2 = 0, is given by C/„ = 14 = <5n,no- In the top left panel of Fig. [T]we present numerical 
results for the stability of this mode in the (u, /o) plane, for a particular value of the lattice coupling constants, 
Ki = K2 = 0.2. Represented in colors is the largest imaginary part of the eigenvalues, the stability region being left 
black. 
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FIG. 1: (Color online) The top top left panel shows the (in)stability regions for the onsite discrete solitons in the (a, /o) plane 
with Ki — K2 ~ 0.2. The top right panel depicts a generic example of the soliton's profile and its stability spectrum at /i = 5.4. 
The bottom panels shows the evolution of the two components \unf and |unp, at n = 25,26 (left) and n = 25 (right), for 
/i = 5.4 and 5.6, as found from simulations of the underlying equations ([T]), which indicates a perfect agreement with the 
threshold value predicted by the averaged equation. 



As seen in the top left panel, for /o = 0.25, the stability window of the localized mode is —0.1072 < a < 0.1084. 
Using Eq. we conclude that the values of /i supporting the stable onsite mode are given by 5.5 < /i < 6.54 for 

= 10. In the top right panel we present the profile of the discrete soliton at /i — 5.4 and its stability spectrum. The 
solution is unstable due to a collision of a pair of eigenvalues at the origin, resulting in a doublet of purely imaginary 
eigenvalues. 

To check the accuracy of the predictions provided by the averaged equations (O, we use the static solution obtained 
from these equations as initial conditions for simulations of underlying equations ([1]), inverting transformation ^ for 
this purpose. Shown in the left bottom panel of Fig. [T] is the generated evolution of the discrete soliton, just near 
the stability border shown in the top right panel of the same figure. One can see that the unstable soliton tends to 
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rearrange itself into a different time-periodic localized solution. We will see later that in the averaged equation ^ 
this corresponds to an asymmetric mode. On the other hand, in the right bottom panel we show the evolution of 
a stable onsite soliton at /i = 5.6, from which one can conclude that our averaged equations accurately predict the 
stability threshold, in comparison with Eqs. ([T]). 



C. Two-site modes 



Proceeding from the fundamental onsite solitons (single-site modes) to the consideration of multi-site ones, a natural 
object is an in-phase two-site state, i.e., an intersite discrete soliton. However, our numerical analysis (not shown here) 
reveals that this mode is unstable, in terms of the averaged equations, for all parameter values. Explicit simulations of 
the evolution of the counterpart of this state within the framework of the original equations ([1]) confirm its instability. 

Next, we consider an out-of-phase fundamental two-site state, alias a twisted localized state, which, in the anti- 
continuum limit, is seeded by ansatz Un = Vn = 5n,n„ — (^n.no+i- Results of the numerical analysis of this family of 
odd discrete solitons are presented in Fig. [2] With ki = K2 = 0.05 and /i = 4, which corresponds to cr = 0.455, the 
mode is unstable in the range of 0.223 < fo < 0.469. In the left and right bottom panels of the figure, we show the 
evolution of the twisted mode for fo — 0.22 and /o — 0.23, respectively. Again we observe a good agreement between 
the stability threshold predicted by averaged equations ^ and the original system, Eqs. ([T]). 
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FIG. 2: (Color online) The top left panel is the same as in Fig. \T\ but for a twisted mode with ki — Hi2 ~ 0.05. The top right 
panel shows the profile of the mode and its stability spectrum for /i = 4 and fo = 0.23. Bottom panels show the evolution of 
stable and an unstable modes at fo = 0.22 and 0.23, respectively. 



IV. SYMMETRY BREAKING: PITCHFORK AND HOPF BIFURCATIONS 



As reported in Refs. [12,] and [ISJ, linearly coupled two-component systems, which are similar to the present one, 
give rise to a symmetry-breaking bifurcation, which generates asymmetric solitons, while destabilizing symmetric 
fundamental ones. The analysis of the present model reveals that a similar bifurcation occurs in averaged equations 
([3]), and it generates the respective asymmetric mode. Presented in the top left panel of Fig. [3] is the (in)stability 
region of the asymmetric mode for ki — K2 — 0.2. In this panel, the white curve designates the bifurcation line, below 
which the asymmetric mode exists. The analysis demonstrates that the bifurcation in the averaged system ^ is of 
the pitchfork type, like in the recently studied related models [H, [l3|)- As a particular example, we depict in the 
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top right panel of Fig. [3] the solution and the stability spectral plane for an unstable asymmetric onsite mode with 
/o — 0.23 and /i = 7, which corresponds to cr = —0.19. 
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FIG. 3: (Color online) The same as in Fig. [T] but for the asymmetric mode, generated by the symmetry-breaking bifurcation, 
which happens along the while curve in the top left panel. See the text for further details. 

It is worthy to note that the stability border of the asymmetric mode in the top left panel exactly corresponds to 
tr = 0, the asymmetric discrete soliton being stable at ct > 0. This observation may be explained by the fact that the 
self-phase-modulation nonlinear terms in Eqs. ([3|) dominate over their cross-phase-modulation counterparts just at 
cr > 0. 

Before we proceed to checking the accuracy of the averaged equations, it is necessary to make the following notice 
concerning transformation ([2]): if and are static, then |u„(i)|^ and |w„(i)p will also be static, provided 

that Uri and Vn are symmetric. On the other hand, if ?7„ and are asymmetric, the transformation yields |u„(t)p and 
|w„(t)|^ oscillating periodically in time. Therefore, the pitchfork bifurcation in the averaged equations ([5]) corresponds 
to a Hopf bifurcation in the original system ([T]). 

Taking the above-mentioned fact into regard, an unstable asymmetric mode in the averaged equations is expected 
to correspond to an unstable periodic solution of Eqs. ([T]). In the bottom left panel of Fig. [31 we display the dynamics 
of the mode shown in the top right panel of the same figure, represented by M25 and W25 . It is seen that the oscillating 
solution is indeed unstable, as predicted by the averaged system. 

As another test of the accuracy of the approximation based on the averaging method, in the bottom right panel of 
Fig. [3] we display the dynamics of a stable asymmetric mode with /i = 6 (cr = 4.15 x 10~^), as produced by direct 
simulations of the original equations ([T]). Again, we observe a reasonable agreement between the two systems. 

V. CONCLUSION 

In this work, we have studied discrete solitons in the system of two linearly coupled DNLSEs (discrete nonlinear 
Schrodinger equations), which describes several physical situations in BEG. For the case of the coupling constant 
rapidly varying in time, we have derive the averaged system of coupled generalized DNLSEs. In addition to the linear 
coupling, the latter system contains nonlinear interaction terms of the four-wave-mixing type. We have constructed 
families of symmetric single-site localized states (fundamental discrete solitons), as well as in-phase and twisted two-site 
modes, and identified their stability regions. In particular, the two-site in-phase states are completely unstable, while 
their twisted counterparts have a stability area. The symmetry-breaking bifurcation of the fundamental discrete 
soliton, and the resulting asymmetric fundamental localized mode have been found too, and the stability of the 
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asymmetric mode was explored. By means of direct simulations, it was verified that, in all the cases, the averaged 
equations provide for a very accurate approximation, in comparison with the underlying system featuring the rapidly 
varying linear coupling. 
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